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Abstract. - We investigate the burst dynamics during drainage going from low to high in- 
jection rate at various fluid viscosities. The bursts are identified as pressure drops in the 
pressure signal across the system. We find that the statistical distribution of pressure drops 
scales according to other systems exhibiting self-organized criticality. The pressure signal was 
calculated by a network model that properly simulates drainage displacements. We compare 
our results with corresponding experiments. 



Since the early 1980s physicists have paid attention to the complex phenomena observed 
when one fluid displaces another fluid in porous media. The papers that have appeared in 
the literature mostly refer to the rich variety of displacement structures that is observed due 
to different fluid properties like flow rate, viscosity, interfacial tension, and wettability. The 
major displacement structures have been found to resemble structures generated by geometrical 
models like invasion percolation (IP) |[ |, |], DLA [g, |, |, §, and anti-DLA || j|. Only a 
few authors |)[ [n| |ll| have addressed the interplay between the displacement structures and 
the evolution of the fluid pressure. In slow drainage when non-wetting fluid displaces slowly 
wetting fluid in porous media, the pressure evolves according to Haines jumps |)| [l2| The 
displacement is controlled solely by the pressure difference between the two fluids across a 
meniscus (the capillary pressure), and the non- wetting fluid invades the porous medium in a 
series of bursts accompanied by sudden negative pressure drops. 

The purpose of this paper is to study the dynamics of the fluid pressure during drainage 
going from low to high displacement rates. To do so, we examine the statistical properties of the 
sudden negative pressure drops due to the bursts. We find that for a wide range of displacement 
rates and fluid viscosities, the pressure drops act in analogy to theoretical predictions of 
systems exhibiting self-organized criticality like IP. Even at high injection rates, where 
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Fig. 1. - The pressure signal as function of injection time, P(t), for one simulation at low displacement 
rate in a narrow time interval. The horizontal line defines the pressure valley of a burst that last a 
time period T = £2 — £i- Note that the valley may contain a hierarchical structure of smaller valleys 
inside. The vertical line indicates the size of a local pressure jump Ap inside the valley. 



the connection between the displacement process and IP is more open, the pressure drops 
behave similar to the case of extreme low injection rate, where IP apply. The pressures are 
calculated by a network model that properly simulates the fluid- fluid displacement. Moreover, 
we measure the fluid pressure in drainage experiments and compare that with our simulation 
results. 

In the simulations a burst starts where the pressure drops suddenly and stops where the 
pressure has raised to a value above the pressure that initiated the burst (see fig. |f|). Thus, 
a burst may consist of a large pressure valley containing a hierarchical structure of smaller 
pressure jumps (i.e. bursts) inside. A pressure jump, indicated as Ap in fig. |], is the pressure 
difference from the point where the pressure starts decreasing minus the pressure where it 
stops decreasing. We define the size of the pressure valley (valley size) to be x = Si^Pi, 
where the summation index i runs over all the pressure jumps Api inside the valley. The 
definition is motivated by experimental work in ref. For slow displacements we have 

that x is proportional to the geometric burst size s, being invaded during the pressure valley. 
This statement has been justified in ref. |l3j, where it was observed that in stable periods, the 
pressure increased linearly as function of the volume being injected into the system. Later, in an 
unstable period where the pressure drops abruptly due to a burst, this volume is proportional 
to s. At fast displacements the pressure may no longer be a linear function of the volume 
injected into the system. Therefore, a better estimate of s there, is to compute the time 
period T of the pressure valley (fig. [I]). Since the displacements are performed with constant 
rate, it is reasonable to assume that T is always proportional to the volume being injected 
during the valley and hence, T cx s. 

We have computed the distributions of \ an d T from the pressure signals of simulations 
and experiments. We find that the distributions are consistent with a power law, independent 
of injection rate and fluid viscosities (figs. ^| and |J) and that the distribution of pressure jumps 
Api, follows an exponential decreasing function (fig. ^). 

The network model used in the simulations is thoroughly discussed in refs. |f5| and ]l6| and 
only its main features are presented below. The porous medium consists of a two-dimensional 
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(2D) square lattice of cylindrical tubes oriented at 45° relative to one of the edges of the lattice. 
Four tubes meet at each intersection where we put a node having no volume. The disorder 
is introduced by moving the intersections a randomly chosen distance away from their initial 
positions, giving a distorted square lattice. The distances are chosen in the interval between 
zero and less than one half of the grid size to avoid overlapping intersections in the new lattice. 
We let dij denote the length of the tube between node (intersection) i and j in the lattice and 
Tij = dij /2a defines the corresponding radius of the tube. Here a is the aspect ratio between 
the tube length and its radius. 

The tubes are initially filled with a wetting fluid of viscosity /i w , and a non- wetting fluid of 
viscosity ^ nw is injected at constant injection rate Q along the bottom row. The wetting fluid 
is displaced and flows out along the top row and there are periodic boundary conditions in the 
horizontal direction. The fluids are assumed incompressible and immiscible and an interface 
(meniscus) is located where the fluids meet in the tubes. The capillary pressures of the menisci 
behave as if the tubes where hourglass shaped with effective radii following a smooth function. 
Thus, we let the capillary pressure p c be a function of the meniscus' position in the tube in 
the following way: p c — (2j/r)[l — cos(2nx/d)]. Here we have omitted the subscripts ij. The 
first term results from Young-Laplace law when assuming that the principal radii of curvature 
of the meniscus are equal to the radius of the tube, and that the wetting fluid perfectly wets 
the medium. 7 denotes the interfacial tension between the fluids. In the second term x is the 
position of the meniscus in the tube, i.e. < x < d. The advantage of the above approach 
is that we include the effect of local readjustments of the menisci on pore level [^5), which is 
important for the description of the burst dynamics ^ [H| . 

The fluid flow c/y through a tube from node i to node j, is solved by using Hagen-Poiseuille 
flow in cylindrical tubes and Washburn's approximation jlTj for menisci under motion giving, 
qij — —(oijkij j '^,ij){pj — pi — Pc.ij) I dij- Here pi and pj are the pressures at the nodes, p c ^j 
is the capillary pressure if one or two menisci are present in the tube, and /iy is the effective 
viscosity of the fluids occupying the tube, fey and cry is the permeability and the cross section 
of the tube, respectively. By inserting the above equation into Kirchhoff equations at every 
node, Y]j q^ — 0, constitutes a set of linear equations which are solved for the nodal pressures 



Pi. The set of linear equations is solved by the Conjugate Gradient method |lq|. See refs. |15 



and |16[ for how the menisci are updated and other numerical details about the network model. 

To characterize the fluid properties used in the simulations, we use the capillary number 
C a and the viscosity ratio M. C a , denoting the ratio of capillary and viscous forces, is in the 
following defined as C a = Q/i/£7. Here [i is maximum viscosity of p nw and /j, w , and £ is the 
cross section of the inlet. The viscosity ratio M, is defined as M = /i nw /fi w . 

We have performed three different series of simulations with M = 0.01, 1, and 100, 
respectively. In each series C a was varied by adjusting the injection rate Q. To obtain reliable 
average quantities we did 10 to 20 simulations of different distorted lattices, at each C a . The 
lattice size of the networks was 60 x 90 nodes for M = 0.01, 40 x 60 nodes for M — 1, and 
25 x 35 nodes for M = 100. In all simulations we set 7 = 30 dyn/cm, and the radii of the 
tubes were inside the interval [0.08, 0.72] mm. The average tube length was always 1 mm. The 
parameters were chosen to be close to the experimental setup in Q . 

For all simulations we calculated the hierarchical valley size distribution iV a n(x). The 
distribution was calculated by including all valley sizes and the hierarchical smaller ones within 
a large valley (see fig. [l]). The result for high, intermediate, and low C a when M — 1 and 
AI = 100 is shown in a logarithmic plot in fig. |[ Identical results were obtained for M = 0.01. 
In order to calculate the valley sizes at large C a , we subtract the average drift in the pressure 
signal due to viscous forces such that the pressure becomes a function that fluctuates around 
some mean pressure. 
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Fig. 2. - The hierarchical valley size distribution iV a ii(x), for simulations between low and high C a 
with M = 1 (O ,□,<>) and M = 100 (A,< , V). The slope of the solid line is -1.9. Inset: The 
cumulative valley size distribution N(x>X*)> f° r bursts that start in a narrow pressure strip for the 
simulation performed at C a = 1.6 x 10 -5 . The slope of the solid line is —0.5. 



By assuming a power law iV a n(x) oc x~ Ta11 our best estimate from fig. || is r a n = 1.9 ± 0.1, 
indicated by the slope of the solid line. At low \ m n g- H typically only one tube is invaded 
during the valley and we do not expect the power law to be valid. Similar results were obtained 
when calculating the hierarchical distribution of the time periods T of the valleys, denoted as 
N aU (T). 

In IP the distribution of burst sizes N(s), where s denotes the burst size, is found to obey 
the scaling relation J9|, [H| [ll| ||(| 

N( S )KS- T 'g(s°(f -f c )). (1) 

Here / c is the percolation threshold of the system and g(x) is some scaling function, which 
decays exponentially when x 3> 1 and is a constant when x — > 0. r' is related to percolation 
exponents like t' = 1 + Df/D — \j(Dv) pQ] , where Df and D is the fractal dimension of the 
front and the mass of the percolation cluster, respectively. Df depends on the definition of 
the front, that is, D{ equals D c for external perimeter growth zone |^l|, ^2) and D\ Y for hull 
perimeter growth zone |2^, |23|. v is the correlation length exponent in percolation theory 
and cr = 1/{vD) j22j. In cq. (|l]) a burst is defined as the connected structure of sites that is 
invaded following one root site of random number /o, along the invasion front. All sites in the 
burst have random numbers smaller than /o, and the burst stops when / > /q, is the random 
number of the next site to be invaded Ell . 



By integrating eq. (|l|) over all /q in the interval [0, / c ] Maslov 1 14 deduced a scaling relation 
for the hierarchical burst size distribution N a \\{s) following 

AT all (s) cx S - T ^, (2) 

where r a u = 2. 

In the low C a regime in fig. ^, the displacements are in the capillary dominated regime and 
the invading fluid generates a growing cluster similar to IP m M, ^5|, [2(| . In this regime we also 
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Fig. 3. - The cumulative pressure jump distribution function N(P > P*), for simulations performed 
with viscosity matched fluids (M = 1). The dashed lines are fitted exponential functions. 



have that x oc s jl3| and hence -/V a n(x) corresponds to N^u(s) in eq. (|^). Thus, in the low C a 
regime we expect that -ZV aU (x) follows a power law with exponent r a n = 2 which is confirmed 



by our numerical results. Similar results were obtained in ref. |13 



The evidence in fig. |2|, that r a n does not seem to depend on C a , is very interesting and 
new. At high C a when M = 0.01 an unstable viscous fingering structure generates and when 
M > 1 a stable front develops. It is an open question how these displacement processes map 
to the proposed scaling in eq. (||) . We note that in the high C a regime the relation x k s may 
not be correct and T is preferred when computing N^n- However, the simulations show that 
N aU (x) ~ even at high C„. 

In JlJ] it was pointed out that T a \\ is super universal for a broad class of self-organized 
critical models including IP. Our result in fig. ^| indicates that the simulated displacement 
processes might belong to the same super universality class even at high injection rates. 

Maslov Ql also calculated the time-reversed (backward) hierarchical burst size distribu- 
tion and predicted that this distribution should follow a power law with a model-dependent 
exponent rh,. In our case we are dealing with 2D IP with trapping giving t^ u — 1.68. We 
have calculated t^j of our simulations by simply reversing the time axis in the pressure signal 
in fig. [l] and repeating the steps which led to fig. From that we obtain r all = 1.7 ± 0.1 which 
is consistent with the predictions in 

In the inset of fig. || we have plotted the cumulative valley size distribution N(x>x*) f° r 
the simulation at lowest C a = 1.6xl0~ 5 with M — 1. N(x>X*) was calculated for bursts 
that starts at pressures in a narrow strip between 2800 and 3100 dyn/cm 2 where 3100 is the 
maximum pressure during the displacement. From eq. (Q) we have that N(s) oc s~ T for 
bursts that start close to the percolation threshold f c . In our simulations f c corresponds to 
the maximum pressure. It is hard to observe any power law in the inset of fig. ^ however, 
if we assume one, our best estimate is 1 — r' = —0.5 as indicated by the slope of the solid 
line. In |H| simulations and experiments gave 1 — t' = —0.45 ± 0.10. We need larger system 
sizes and more simulations to improve our statistics, but we conclude that our result are in 
agreement of [|l3[ . 

We have also calculated the cumulative pressure jump distribution function N(P> P*) for 
the simulations with M = 1 and 100 at various injection rates. Here P = Ap/ (Ap) where 

). The 
Both 



(Ap) is the mean of the local pressure jumps Ap in the pressure signal (see fig 
result for two simulation, one at high and the other at low C a , is plotted in fig 
were performed with viscosity matched fluids (M = 1). The distributions have been fitted 
to exponentially decreasing functions drawn as dashed lines in fig. [| At low C a we find 
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Fig. 4. - The hierarchical distribution -/Vaii(T) of the valley time T during a burst, for experiments 
(open symbols) and simulations (filled symbols) at various C a with M = 0.017 and M = 0.01, 
respectively. The slope of the solid line is —1.9. 



N(P> P*) oc e~ 138F , which is consistent with results in At high C a the distribution 
function was fitted to e _102P . The pre-factor in the exponent of the exponential function 
seems to change systematically from about 1.4 to 1.0 as C a increases. Similar results were 
obtain from simulations performed with M = 100. 

We have performed four drainage experiments where we used a 110 x 180 mm transparent 
porous model consisting of a mono-layer of randomly placed glass beads of 1 mm, sandwiched 
between two Plexiglas plates || . The model was initially filled with a water-glycerol mixture 
of viscosity 0.17 P. The water-glycerol mixture was withdrawn from one of the short side of 
the system at constant rate by letting air enter the system from the other short side. The 
pressure in the water-glycerol mixture on the withdrawn side was measured with a pressure 
sensor of our own construction. 

From the recorded pressure signal we calculated the hierarchical distribution of time periods 
of the valleys, N a \\(T). At low C a this corresponds to iVail(s) m eq. @). Because of the relative 
long response time of the pressure sensor, rapid and small pressure jumps due to small bursts 
are presumably smeared out by the sensor and the recorded pressure jumps are only reliable 
for larger bursts. Hence, from the recorded pressure signal T appears to be a better estimate 
of the burst sizes than \. 

In fig. [I] we have plotted the logarithm of AT aU (T) for experiments (open symbols) and 
simulations (filled symbols) performed at four different C a , respectively. To collapse the 
data iV a ii(T) and T were normalized by their means. In the simulations M = 0.01 while 
in the experiments M = 0.017 where we have assumed air to have viscosity 0.29 x 10~ 2 P. We 
observe that the experimental result is consistent with our simulations and we conclude that 
N aU (T) oc r 19±01 . This confirms the scaling of AT a u(x) in fig. |. We have also calculated the 
time- reversed distribution of N a \\(T) and the result of that is consistent with the time-reversed 
distribution that was calculated of the simulations in fig. |^. 

Note that when comparing the C a 's of the experiments with the ones of the simulations in 
fig. ^, we have to take into account the different system sizes. The length of the experimental 
setup is about three times larger than the length of the simulation network. Therefore we 
expect that in the experiments, viscous fingering develops at C a 's of about three times less 
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than in the simulations. 

In summary we find that T a n = 1.9 ± 0.1 for all displacement simulations going from low to 
high injection rates when M = 0.01, 1, and 100. This is also confirmed by drainage experiments 
performed at various injection rates with M = 0.017. At low injection rates the result is 
consistent with the prediction in |l4| (rail = 2), which was deduced for a broad spectrum 
of different self-organized critical models including IP. The evidence that r a n is independent 
of the injection rate, may indicate that the displacement process belongs to the same super 
universality class as the self-organized critical models in fl4|| , even where there is no mapping 
between the displacement process and IP. The good correspondence between our simulation 
results and the drainage experiments in fig. || and also the results reported at slow drainage 
in |l^| , demonstrates that the burst dynamics is well described by our simulation model. 

The authors thank S. Roux for valuable comments. The work is supported by the Norwegian 
Research Council (NFR) through a "SUP" program and we acknowledge them for a grant of 
computer time. 



REFERENCES 

Wilkinson D. and Willemsen J. F., J. Phys. A, 16 (1983) 3365. 
Lenormand R. and Zarcone C, Phys. Rev. Lett, 54 (1985) 2226. 
Cieplak M. and Robbins M. O., Phys. Rev. Lett, 60 (1988) 2042. 
Witten T. A. and Sander L. M., Phys. Rev. Lett, 47 (1981) 1400. 
L. Paterson, Phys. Rev. Lett, 52 (1984) 1621. 
Chen J.-D. and Wilkinson D., Phys. Rev. Lett, 55 (1985) 1892. 
Maloy K. J., Feder J. and JOSSANG T., Phys. Rev. Lett, 55 (1985) 26881. 
Lenormand R., Touboul E. and Zarcone C, J. Fluid Mech., 189 (1988) 165. 
Maloy K. J., Furuberg L., Feder J. and Jossang T., Phys. Rev. Lett., 68 (1992) 2161. 
van der Marck S.C., Matsuura T. and Glas J., Phys. Rev. E, 56 (1997) 5675. 
van der Marck S. C. and Glas J., Eur. J. Mech., B/Fluids, 16 (1997) 681. 
W. B. Haines, J. Ayr. Set., 20 (1930) 97. 

Furuberg L., Maloy K. J. and Feder J., Phys. Rev. E, 53 (1996) 966. 
S. Maslov, Phys. Rev. Lett, 74 (1995) 562. 

Aker E., Maloy K. J., Hansen A. and Batrouni G. C, Transp. Porous Media, 32 (1998) 
163. 

Aker E., Maloy K. J. and Hansen A., Phys. Rev. E, 58 (1998) 2217. 
E. W. Washburn, Phys. Rev., 17 (1921) 273. 
Batrouni G.G. and Hansen A., J. Stat. Phys., 52 (1988) 747. 

Sapoval B., ROSSO M., and Gouyet J. F., in Fractals' Physical Origin and Properties, edited 
by L. Pietronero, (Plenum Press, New York) 1989. 

Martys N., Robbins M. O. and Cieplak M., Phys. Rev. B., 44 (1991) 12294. 
Grossman T. and Aharony A., J. Phys. A: Math. Gen., 20 (1987) L1193. 
Stauffer D. and Aharony A., Introduction to percolation theory, (Taylor & Francis, London, 
Great Britain) 1992. 

R. F. VOSS, J. Phys. A: Math. Gen., 17 (1984) L373. 
Roux S. and Guyon E., J. Phys. A: Math. Gen., 22 (1989) 3693. 
de Gennes P.G. and Guyon E., J. Mec. (Pans), 17 (1978) 403. 

Chandler R., Koplik J., Lerman K. and Willemsen J. F., J. Fluid Mech., 119 (1982) 249. 



